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ABSTRACT 

The merging history of dark matter halos is computed with the Merging Cell Model 
proposed by Rodrigues & Thomas (1996). While originally discussed in the case of 
scale-free power spectra, it is developed and tested here in the framework of the cold 
dark matter cosmology. The halo mass function, the mass distribution of progeni- 
tors and child halos, as well as the probability distribution of formation times, have 
been computed and compared to the available analytic predictions. The halo auto- 
correlation function has also been obtained (a first for a semi-analytic merging tree) , 
and tested against analytic formulae. An overall good agreement is found between 
results of the model, and the predictions derived from the Press & Schechter theory 
and its extensions. More severe discrepancies appear when formulae that better de- 
scribe N-body simulations are used for comparison. In many instances, the model can 
be a useful tool for following the hierarchical growth of structures. In particular, it is 
suitable for addressing the issue of the formation and evolution of galaxy clusters, as 
well as the population of Lyman-break galaxies at high redshift, and their clustering 
properties. 
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In the basic picture of the formation of cosmic structures, 
the Universe is dominated by a dark matter (DM) com- 
ponent, and small perturbations in the initial density field 
grow in amplitude proportionally to a linear growth factor, 
until they approach unity. Then, non-linear effects domi- 
nate their evolution, and the regions stop expanding with 
the Universe, collapse, and virialise, thus forming DM halos. 
In hierarchical scenarios, like the cold dark matter (CDM) 
model, small-scale inhomogeneities collapse first, and then 
aggregate via merging to generate larger structures. Since 
galaxies form by the collapse and cooling of baryonic gas 
within DM halos, and their history is greatly influenced by 
that of their surrounding halos (e.g., Lemson & Kauffmann 
1999), it is important to understand how these objects form 
and evolve with time. 

The most realistic way for following the history of DM 
halos is by means of N-body simulations, but they require 
huge amounts of RAM memory and are computationally 
expensive. Therefore, they are often limited to a modest 
dynamic range, and to few different cosmological scenarios. 

The simplest alternative approach is to consider only 
the linear regime of growth of density fluctuations, and de- 



scribe the non-linear evolution and collapse by means of 
the spherical 'Top-Hat' model (Gunn & Gott 1972). In this 
formalism, the formation of a DM halo of mass M at red- 
shift z is described by identifying in the initial density field 
smoothed on a scale M, and linearly extrapolated to red- 
shift z, a region having overdensity equal to a given thresh- 
old value. Starting with Gaussian initial conditions, Press 
& Schechter (1974, hereafter PS) interpreted the probabil- 
ity of finding such regions as the number density of halos 
of mass M, that formed at redshift z (see Sec. 3.1). The 
PS mass function has been extensively tested against N- 
body simulations, and has been found to be in reasonably 
good agreement with numerical results (e.g., Efstathiou et 
al. 1988; Carlberg & Couchman 1989; Lacey & Cole 1994; 
Gelb & Bertschinger 1994). However, systematic deviations 
both at small and high masses have been recognised, with 
the PS formula predicting too many low mass halos, and un- 
derestimating the number of massive objects (e.g., Jain & 
Bertschinger 1994; Gross et al. 1998; Somerville et al. 1998; 
Tormen 1998; Lee & Shandarin 1999; Sheth & Tormen 1999; 
and references therein). A better agreement with numerical 
results is obtained when ellipsoidal rather than spherical col- 
lapse models are considered (Monaco 1995, 1997a, b; Bond 
& Myers 1996; Audit, Teyssier & Alimi 1997, 1998; Lee & 
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Shandarin 1998, 1999, hereafter LS98 and LS99; Sheth & 
Tormen 1999; Sheth, Mo, & Tormen 1999; and references 
therein). 

Extensions of the PS theory (Bower 1991; Bond et al. 
1991) follow the redshift evolution of the halo population as 
a whole, by deriving the conditional probability of finding 
progenitors of mass M p at redshift z p , given their child halos 
of mass M at z , and vice versa (see Sec. 3.2). By means 
of the extended Press & Schechter (EPS, hereafter) theory, 
the distribution of halo formation and survival times, as well 
as their merger rate, can also be obtained (Lacey & Cole 
1993, 1994, hereafter LC93 and LC94; see also Sec. 3.3). 
The few comparisons between these analytic predictions and 
numerical results reveal a general good agreement, even if 
discrepancies similar to those of the mass function have been 
pointed out (LC94; Somerville et al. 1998; Tormen 1998). 

Still based on the EPS theory, analytic predictions for 
halo bias in the Lagrangian space of initial conditions have 
been obtained (Mo & White 1996, MW96 hereafter; Catelan 
et al. 1998, CLMP; Porciani et al. 1998; Sheth & Lemson 
1999b; Sheth, & Tormen 1999; Sheth, Mo, & Tormen 1999; 
and references therein). The halo auto-correlation function 
£,hh(y) is then the product of the halo bias with the corre- 
lation function of the underlying matter. The predicted £hh 
is in good agreement with that in N-body simulations for 
massive objects, but its amplitude is too large for low-mass 
halos (Porciani, Catelan & Lacey 1999, PCL hereafter; Jing 
1999; Sheth & Tormen 1999). However, Jing (1999) propose 
an empirical fitting formula (see Sec. 3.5) that provides a 
good description of halo clustering on the whole range of 
masses (see also Sheth & Tormen 1999; Sheth, Mo, & Tor- 
men 1999). 

While the PS and EPS formalisms describe the mean 
statistical properties of the population as a whole, sev- 
eral models of the individual merger history of DM halos 
have been proposed (Cole & Kaiser 1988; Kauffmann & 
White 1993, KW93 hereafter; Rodrigues & Thomas 1996, 
RT96 hereafter; Somerville & Kolatt 1999; Sheth & Lem- 
son 1999a). Each model presents some advantages and some 
drawbacks with respect to the others. For example, the 
"block model" of Cole & Kaiser (1988) partly retains the 
spatial information, but it is affected by the discretisation 
of both halo masses (in powers of two), and positions. The 
KW93 merging tree presents a more continuous spectrum of 
masses, but a grid of collapse redshifts is imposed, and the 
relative positions of halos are unknown. Moreover, it repro- 
duces exactly the mean progenitor mass distribution, but 
mass conservation is enforced only approximately, while the 
opposite holds in the model of Somerville & Kolatt (1999). 

In this paper, we focus on the "Merging Cell Model" 
(MCM) Q proposed by Rodrigues & Thomas (1996), which 
has the same characteristics of simplicity and speed as the 
other merging tree algorithms, but also presents some ma- 
jor advantages. Since it is based on an actual realisation of 
the initial density field, it is much closer to the spirit of 
N-body simulations, thus allowing direct comparisons with 
numerical results, and it also seems to take into account 
the spatial correlations of density fluctuations (Nagashima 
& Gouda 1997) . Moreover, no specific collapse times are im- 



as first coined by Nagashima & Gouda (1997). 



posed a priori, and halos form with a continuous spectrum 
of masses, and a variety of (Lagrangian) shapes. Also the 
spatial information about the relative location of halos is 
retained by construction, thus allowing to study their clus- 
tering properties. While in the original paper, the authors 
only discuss the halo mass function in the case of scale-free 
power spectrum, here we consider the more realistic stan- 
dard CDM (SCDM) cosmology (see Sec. 3). Moreover, we 
test the model reliability also in terms of the mass distri- 
bution of progenitor and child halos, the behaviour of the 
largest progenitor mass as a function of redshift, and the 
probability distribution of formation times. For the first time 
for a semi-analytic merging tree, also the halo two-point 
correlation function is computed, and we test it against the- 
oretical predictions. We outline the method in Sec. 2, define 
these quantities and compare them to the analytic predic- 
tions in Sec. 3. Discussion and conclusions are presented in 
Sees. 4 and 5, respectively. 



2 THE ALGORITHM 

2.1 Basic Principles 

At an 'initial time' ti, consider the density field p(x,ti) of 
the Universe characterised by a mean value ~p(ti), and small 
perturbations <5(x,i ; ) = p(x,£ ; )/p(ti) — 1. Through grav- 
itational instability, the amplitude of density fluctuations 
start growing proportionally to a linear growth factor D(t), 
i.e., <5(x,£) = S(x,t[) x D(t)/D(ti). Such a linear growth 
law strictly holds only when perturbations are much smaller 
than unity, but it is useful to extrapolate it also into the non- 
linear regime. In fact, the 'Top-Hat' model (Gunn & Gott 
1972) shows that the formation of a bound virialised object 
of mass M occurs at the time U when the density contrast of 
a spherical region in the initial density field, smoothed at a 
scale M, reaches a critical value <5 C . This in turn corresponds 
to a value <5j. ln (£f) of the density field linearly extrapolated 
to that time, or to a value S^ n (to) if the extrapolation is car- 
ried on until the present epoch to- For an Einstein-de Sitter 
Universe, D(t) oc (1 + z)' 1 , 5 C ^ 178, <5j? n (z f ) ~ 1.686, and 



Sr(z = 0) 



i (1 + Zt). It is therefore sufficient to know 



the values of the density field linearly extrapolated to z — 0, 
for determining the formation epochs of DM halos. 

Because this is the approach we adopt in the paper, we 
choose henceforth to change the notation and we denote the 
density field linearly extrapolated to z = as 8. Therefore, 
8 C = 1.686, and the collapse redshift of halos is zi = 5/8 c — 1. 



2.2 The Method 

The MCM is based on an actual realisation of the density 
field, obtained with a standard initial condition generator, 
by Fourier transforming waves of random phase and ampli- 
tude drawn from a Gaussian distribution of zero mean, and 
variance given by the chosen power spectrum. 

A value of the density contrast S is assigned to each of 
the L 3 base cells (bes) composing a periodic cubic box of 
side L (for simplicity, L=2 ! , where I is a positive integer). 
Density fluctuations are then averaged within cubic blocks 
of side 2, 4, 8, ..., L. At each of these smoothing levels, 
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a set of 8 overlapping grids displaced one relative to an- 
other by half a block length in each coordinate direction is 
used. This ensures that the density peaks are always approx- 
imately centred within one of the blocks in each smoothing 
hierarchy. At this point, one has a total of (15 L 3 — 8)/7 base 
cells and cubic blocks, with side ranging from 1 to L, mutu- 
ally overlapping in the volume of the box. Each of them is 
characterised by a value of the density contrast 5. 

All base cells and blocks are then ordered in a single 
list in terms of decreasing 8 (or, correspondingly, decreasing 
collapse redshift). The largest value of S in the list fixes the 
earliest Zi of the realisation. It usually corresponds to a base 
cell, which thus becomes the first collapsed object (i.e., the 
first halo). All the elements of the list are then analysed 
one after the other from early times to the present, and the 
specific base cell or cubic block under investigation is called 
investigating region. Whether the investigating region can 
collapse and give rise to a new halo or not is decided by the 
following rules: 

(i) an investigating region that does not overlap with any 
other pre-existing halo collapses and forms a new halo; 

(ii) if there exist two halos, each of them containing half 
of the investigating region, the latter cannot collapse. This 
is to avoid the formation of very elongated structures in 
linking together adjacent halos without the collapse of any 
new matter. If instead there exists only one halo containing 
half (or more) of the investigating region, the latter collapses 
and merges with it, thus forming a new halo; 

(iii) after taking into account condition (ii) , if the inves- 
tigating region overlaps with at least half of one (or more) 
pre-existing halo(s), it collapses and merges with it (them), 
thus forming a new halo. 

Note that in the Lagrangian space of initial conditions, 
mass and volume are equivalent quantities (M=p V), and 
'merging' together the investigating region with one or more 
pre-existing halos does not mean summing up their masses. 
Instead, the mass (volume) of the new resulting halo is that 
of the old ones plus the fraction of the investigating region 
that does not overlap with any already pre-existing object. 

Thanks to the use of overlapping grids and merging cri- 
teria, halos of a large variety of shapes and masses are ob- 
tained. The model also contains information on the relative 
locations of halos, since their positions within the box are 
known, and the effects of discretisation are expected to be 
smaller than in the block model. Moreover, because halo 
formation times are given by the density contrasts in the 
list, they span a continuous range of values. As a drawback, 
the overdensity of the investigating region (which fixes the 
collapse redshift of the new forming object) is not necessar- 
ily equal to the mean overdensity (5h) of the resulting halo. 
Therefore, the position of an object in the merging tree oc- 
casionally differs from that predicted by the linear theory, 
i.e., the assigned zt is not exactly equal to Sh/S c — 1, as it 
should (see Sec. 3.2 in RT96). Finally, the 'linking' and the 
'overlapping' conditions [criteria (ii) and (iii), respectively] 
are reasonable but arbitrary, and different choices would re- 
sult in different mass functions, as discussed by Nagashima 
& Gouda (1997). Yet, it is not clear which are the 'best' 
criteria. Thus, we will adopt the original conditions (ii) and 
(iii) throughout the paper. 



3 TESTS OF THE ALGORITHM 

The MCM is based on the linear theory of growth of den- 
sity fluctuations, and it uses simplified criteria to describe 
the formation and merging history of DM halos. It is there- 
fore necessary to test its reliability by comparing its results 
against those of N-body simulations that directly take into 
account the gravitational interactions between DM particles, 
and are much more realistic in following the dynamics in the 
non-linear regime. 

The model is required to correctly describe not only 
the population of halos at a given redshift, but also how 
this population evolves with time. 

As a first test, the cumulative and the differential mass 
functions in the case of a scale-free power spectrum with 
spectral index n = 0, and —2 have been computed and com- 
pared to those in the original paper (Figs. 4 and 5 in RT96). 
A remarkable agreement has been found. 

Here we consider the SCDM cosmology, and perform 
several tests against the available analytic formulae, to ver- 
ify the reliability of the model results. We set the Hubble 
constant to H = 100 ftkm s -1 Mpc -1 , h = 0.5. The total 
and baryonic density parameters are f2o = 1 and fib = 0.05 
respectively, while that corresponding to the cosmological 
constant is £7a = 0. The transfer function of Bardeen at al. 
(1986) is adopted, and the power spectrum is normalised 
so that the mass variance on scale 8 /i _1 Mpc is equal to 
a 8 = 0.67. 

For a 256 3 base cell realisation in a cubic box of L=100 
Mpc side (i.e., 50/i _1 Mpc, and a total mass of about 3.5 x 
10 16 h~ 1 M Q ), the base cell mass is of about 2x 10 9 ft _1 M , 
and the resulting most massive halo typically has a mass of 
about 5x 10 14 r'M Q (but see Sec. 3.1). The CPU time on 
a 500 MHz DEC Alpha workstation is only roughly 25 min, 
and typically 700 MB of RAM memory are required. 

In the following sections, results averaged over 10 dif- 
ferent realisations are presented, and error bars correspond 
to the standard deviations of the 10 run sample. 



3.1 Mass Function 

The differential mass function of halos, is defined as the 
comoving number density of halos with mass in the range 
[M, M+dM] at redshift z. This is shown for logarithmic mass 
interval in the histrograms of Fig. hi (left-hand panels), for 
z = and z — 3. 

Also shown for comparison as dotted lines are the cor- 
responding predictions of the PS theory (e.g., LC94): 
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where S c (z) — 8 C x (1 + z), and <r(M) is the mass variance of 
the linearly extrapolated (to z — 0) density field smoothed 
on scale M. In this paper, er(M) is always computed by a fit- 
ting formula analogous to that proposed by White & Frenk 
(1991), with errors smaller than 8% on mass scales ranging 
from 1O 9 M to 1O 15 M . 

Right panels show the cumulative mass fraction for the 
same redshifts, i.e. the fraction of the total mass which is in 
halos of mass above M, at redshift z. 

At z — the overall agreement is good over a large 
range of masses, but a lack of objects at the two ends of the 
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mass function is evident. At small masses, the problem seems 
to be inherent to the method, since this is also the case for 
scale-free power spectra (RT96; Nagashima & Gouda 1997). 
This is a drawback of the model, which limits the reliable dy- 
namical range, and it is probably due to the adopted criteria 
for the formation and merging of halos. In fact, as detailed 
in Sec. 2, the "bricks" for the construction of halos are base 
cells and blocks composed by 8* bcs (i = 1,2, ..). Thus, an 
object of less than 8 bcs can only result if an investigating 
region partly overlaps with a pre-existing halo, but does not 
merge with it, and the fraction of its non-overlapping vol- 
ume (which gives rise to the new halo) is less than 8 bcs. In 
practice, the model requires that some "particular" condi- 
tions happen in order to form halos of mass between 2 and 7 
bcs, thus explaining the underproduction of this kind of ob- 
jects in the resulting mass function. The lack of high-mass 
structures instead, is partly inherent to the method, partly 
due to a statistical fluctuation (in different realisations in 
fact the problem is more or less severe). At high redshift the 
model always tends to produce a larger number of interme- 
diate mass halos, and less massive objects than predicted by 
the PS theory. 

These discrepancies appear to be even more severe if 
compared to results of N-body simulations. It has been 
recently shown, that the PS mass function already tends 
to predict fewer high-mass halos, and more low mass ob- 
jects than those found in the simulations (e.g., Jain & 
Bertschinger 1994; Gross et al. 1998.; Somerville at al. 1998; 
Tormen 1998; LS99; Sheth & Tormen 1999). An analytic 
formula which better agrees with numerical results has been 
obtained by LS98, based on a nonspherical model for the 
collapse of a perturbation, in the frame of the Zel'dovich 
approximation (but see also Monaco 1995, 1997a,b; Audit, 
Teyssier & Alimi 1997; Bond & Myers 1996; Sheth & Tor- 
men 1999; Sheth, Mo, & Tormen 1999). 

In this formalism, the displacement of a particle due 
to the surrounding density field, is simply computed from 
the perturbation potential \& generated by the distribution 
of particles in the initial conditions. The mass density can 
therefore be expressed as a function of the three eigenvalues 
of the deformation tensor (defined as the second derivative 
of ^J), and a virialised bound object forms when the smaller 
one (A3) is positive. The idea is therefore to substitute the 
collapse condition of the spherical Top-Hat model (5 — 5 C ), 
with an analogous one for A3: a DM halo of mass M forms 
when this eigenvalue reaches a critical value A3 C in a region of 
the linearly extrapolated density field, smoothed on a scale 
M. The resulting mass function is: 



dn LS ,„ v _ 25V10 p 
dlnM 1 ,Z} - 2^F M 

with x — A3c(z)/cr, and: 
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where a — <r(M), erfc(x) is the complementary error func- 



Figure 1. Left panels: comoving number density in units of 
Mpc — 3 , per logarithmic mass interval, of halos of mass M at red- 
shifts z = (upper panel) and z = 3 (lower panel), as a function 
of M/Mq. Masses are also shown in units of base cells on the top 
of the figure. Results averaged over 10 realisations of the MCM 
are displayed in the histograms. Error bars show the standard 
deviations of the 10 run sample. The Press & Schechter (1974) 
predictions [eq. (1)] and the Lee & Shandarin (1998) mass func- 
tion [eq.(2)] are plotted as dotted and solid curves, respectively. 
Right panels: cumulative mass fraction of halos with mass larger 
than M at redshifts z = and z = 3, as a function of M. Results 
from the MCM (thick solid line) are compared to the PS and the 
LS98 predictions (dotted and solid curves, respectively). 



tion, and the critical value for A3 has been empirically chosen 
to be A 3c (z) =0.37(1 + 2). 

Figure hi shows that in comparison to the LS98 mass 
function (solid curves) , the MCM presents an excess of small 
objects, and a significant underproduction of high-mass ha- 
los, especially at high redshift. 



3.2 Conditional Mass Functions 

In this section we analyse how the population of halos iden- 
tified at a given time has changed with respect to a different 
epoch. 

Figure @ shows the mass fraction of halos of mass M at 
redshift z , that has already settled at redshift z p in progen- 
itors with masses between M p and M p + dM p . Child halos 
have been selected at z a — and have masses M in the 
range M o /M = [10\ 10 i+1 [, where, from the top to the 
bottom panel, % = 11, 12, 13, 14. The mass distribution of 
their progenitors at Z p — 1 is shown in the left panels, that 
for 2 P = 3 is plotted on the right-hand panels. 

The analytic prediction for the distribution of progeni- 
tor masses is given by (e.g., Bower 1991, Bond et al. 1991): 
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(<5 cp — 5 co ) 



2 (^ -erg) 



(3) 



where <5 c j = 5 C (;£.,•), ctj = a(Mj). 

Because numerical results shown in Fig. H have been 
obtained for ranges of M OJ the analytic formula has been 
computed in two different ways for an accurate comparison. 
In the MCM, the mass of a halo corresponds to the number 
of base cells that compose it. Thus, when expressed in these 
units, the mass M can only assume TV integer values in the 
range [M;,M S ], with TV = (M s - M; + 1): i.e., M = {M ofc = 
Mi + k - 1, k = 1, TV}. Given the number N(M ok ) of child 
halos with mass equal to M fc for each of its possible TV 
values, the mean weighted mass in the range is: 



M = 



fc=i 

N 



M ofc 



(4) 



where N is the total number of halos in the chosen range of 
Mo. Dotted curves in Fig. pj show results from eq.(3) com- 
puted for M = M (whose numerical values are listed in 
the figure caption). The abrupt fall down of dotted curves 
occurs at values of M p near M , because the progenitor mass 
obviously cannot be larger than that of its child halo. Since 
M is lower than M a in each panel, this explains why dotted 
curves are not as extended in M p as the histograms are. 

For a better comparison between numerical and ana- 
lytic results, we have also computed the average progenitor 
mass distribution by summing up the TV single-mass distri- 
butions, each weighted with the fraction of mass in halos of 
mass M fe [N(M ok ) x M ofc ]: 



hf 



d/ 



din M t 



J2 f(M ok )N(M ok )M ok 

fc=l 

J2 TV(Mo fc ) M ofc 



(5) 



where /(M D fc) = (d//dlnM p ) is the single-mass progenitor 
distribution for child halos of mass M fc as given in eq.(3). 
M p is obviously also a function of M , and it is required that 
M p < M . The average progenitor mass distribution is plot- 
ted in Fig. Has solid curves. By construction, (d//dlnM p ) 
is a sum of curves of the same kind of the dotted lines, with 
the sharp cut-off at progenitor masses very similar or equal 
to that of their child halos. This is the reason for the oscil- 
lations in the solid curves for values of M p between Mi and 
M s . 

A lack of objects with masses between 2 and about 
7 bes is apparent, just as was the case for the mass func- 
tion. The MCM also appears to systematically underpro- 
duce progenitors with mass similar to that of their child 
halos. At intermediate masses, an overall good agreement 
between the MCM results and the analytic predictions is 
found, with a possible slight overproduction of halos in the 
model. When compared to N-body simulations, these dis- 
crepancies may become more severe, since simulations ap- 
pear to have fewer/more halos than predicted by EPS theory 
in the intermediate/high mass range (Somerville et al. 1998; 
Tormen 1998). For the less massive child halos (top panels), 
numerical and analytic results only agree over a small range 
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Figure 2. Progenitor mass distribution at redshifts z p = 1 {left 
panels), and z p = 3 (right panels), for child halos in four differ- 
ent mass ranges M at z = (from top to bottom: Mo/Mq = 
[W i ,W i+1 [,i = 11,12,13,14). Average results from 10 realisa- 
tions of the MCM are displayed in the histograms, and their 
standard deviation is shown as error bars. The average number 
of child halos found in each mass range is, from top to bottom: 
N = 34839, 5545, 822, and 67. Dotted curves correspond to the 
progenitor distribution computed from eq. (3) for M equal to the 
mean mass in the corresponding range (M = M ~ 2.7 X 10 , 
2.5 Xl0 12 , 2.3 Xl0 13 , and 1.5xlO 14 M ). Solid curves are the 
average distribution in each mass range, computed from eq.(5). 



of M p . This is due to the lack of low mass objects, and to 
the fact that the minimum mass in the model is limited to 
1 base cell, thus not allowing to accurately follow back in 
time the past history of small halos. 

The reverse conditional probability that a halo of mass 
M p at z p is incorporated at a later time z in a halo of mass 
between M and M + dM , is shown per logarithmic mass 
interval in the histograms of Fig. til Progenitors of mass 
M p /M = [10*, 10 I+1 [, with i = 11,12,13, are selected at 
z p — 1 (left panels) and z p = 3 (right panels), and the mass 
distribution is computed for their child halos at redshift zero. 
No results for progenitors with masses between 10 14 and 
1O 15 M0 are shown, because too few of them have already 
formed at redshift 1 and 3. 

Given all the objects of mass M p at z p , the analytic 
prediction from the EPS theory for the mass distribution of 
their child halos at z , when expressed per mass logarithmic 
interval, is given by (e.g., LC94): 



d/ nu iA/r n M ° Sc ° ( 5c p 
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exp 
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2ct 2 ct 2 (<t 2 



(6) 



where the notation is the same as in eq.(3). 

As before, eq.(6) has been computed in two different 
ways, in order to get an accurate comparison with the 



numerical distributions. Results for M p = 



M p , the mean 



weighted progenitor mass in the range [Mj,M s ] (analogous 
to Mo), are plotted in Fig. ta as dotted lines. A sharp cut 
off occurs for values of M near to M p , because child ha- 
los cannot be less massive than their progenitors. Since M p 
is larger than Mi in each panel, this explains why dashed 
lines are not as extended in M as histograms are. A more 
appropriate conrparison between numerical and theoretical 
results is obtained if the average child mass distribution 
(d//dlnM Q ) is considered, instead of that relative to pro- 
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Figure 3. Mass distribution of child halos at redshifts z = 0, 
given the progenitors at z p = 1 (left panels), and z p = 3 (right 
panels) with masses, M p /M = [10 i ,l0 i+1 [,i = 11,12,13, from 
top to bottom. Average results from 10 realisations of the MCM 
are plotted as histograms, and their standard deviation is shown 
as error bars. The number of progenitors found in each mass range 
at z p = 1 is, from top to bottom: N p = 59898, 6696, and 431, 
while at z p = 3, N p = 52213, 1724, and 8. Dotted curves refer to 
the children mass distribution computed from eq.(6) for the mean 
progenitor mass in the corresponding range (from top to bottom): 
M p ~ 2.6xlO n ,2.2xl0 12 , and 1.7xl0 13 M Q for progenitors at 
z p = 1; M p ~ 2.3xl0 11 ,1.9xl0 12 , and 1.5xl0 13 M Q in the right 
panels. Solid curves are the average distribution in each mass 
range, computed in the same way as in eq. (5), as detailed in the 
text. 



genitors with mean mass M p : d/(M p )/dln(M ). The com- 
putation of (d//dlnM ) is analogous to that in eq.(5), and 
results are plotted in Fig. pi as solid lines. 

An overall agreement is found between MCM results 
and the EPS theory predictions, that in turn fit reasonably 
well N-body simulations (LC94). However an oscillating be- 
haviour of the child halos mass distribution can be recog- 
nised in the histograms. Actually it is more evident when a 
different binning is used (here results are binned on a mass 
grid M = 2 z bcs, i = 0,1,2..), and it seems inherent to 
the method. Also the halo mass function and the progenitor 
distribution present analogous features, and oscillations ap- 
pear to occur with peaks corresponding to the block masses 
of 8 1 , i — 1,2,.. bes, and with troughs in between. Moreover, 
the same trend is found in the halo mass function for the 
scale free power spectrum (in particular for the spectral in- 
dex n = —2, that fits the CDM spectrum over a significant 
range of masses; see RT96). 



3.3 Largest progenitor history 

By analysing the variation with redshift of the largest pro- 
genitor mass, information can be obtained on how halos 
build up in time, whether they preferably form via a contin- 
uous and slow accretion of small objects, or whether their 
mass suddenly increases because of nearly equal-mass merg- 
ing events, or by mergers of several sub-units at the same 
time. A different behaviour is expected for halos of differ- 
ent masses, with larger objects preferably assembling at re- 
cent epochs, and smaller halos showing a more delayed and 
smooth evolution with time. This is shown for instance, in 
KW93, both from their merging tree model, and from N— 
body simulations (their Figs. 5 and 6, respectively). 

We have looked at the past history of halos with cur- 



Figure 4. History of the most massive progenitor of 30 halos 
selected at z = 0. The ;/-axis represents the ratio of progenitor 
mass M p to final halo mass M . The four panels refer to four 
different values of M o /M : 10 11 , 10 12 , 10 13 , 10 14 . 



rent mass M = 10 11 , 10 12 , 10 13 , 1O 14 M , randomly selecting 
30 objects for each value of M to show the scatter in the 
merging histories. The ratio between the mass of the largest 
progenitor Mi, and that of its child halo M is plotted in 
Fig. [4 as a function of 1 + 2. For all the masses, the expected 
trends are obtained, with larger halos preferably assembling 
through major mergers at low z, and smaller objects grad- 
ually forming in a smoother way by accreting small mass 
objects over a larger interval of time. A qualitative good 
agreement of both trends and scatters is also found between 
the present results and those of KW93. Moreover, halo col- 
lapse occurs at more recent epochs here, as expected when 
a SCDM cosmology is considered instead of an open model 
(fio = 0.2 in KW93). 



3.4 Formation redshift 

In the hierarchical clustering scenario, massive halos form 
by accretion of lower mass structures. Therefore, their for- 
mation redshift is expected on average to be lower than that 
of small objects. Actually, because of the continuous evolu- 
tion in mass due to the hierarchical nature of the process, 
the definition of 'halo formation time' is not straightforward. 
In this paper, we adopt the definition of LC93, as the time 
when half the mass of the halo is assembled, i.e., when a 
progenitor with mass equal to half or more that of its child 
halo appears for the first time. 

Figure H shows the distribution of formation redshift 
for halos with mass M = [10*, 5x10*] M , i = 11, 12, 13, 14 
at z — 0. In agreement with results of the previous section, 
high mass objects tend to form at more recent epochs, while 
lower mass halos typically collapse earlier and over a larger 
interval of time. The mean formation redshifts for halos in 
the four mass ranges, from lower to higher M , are: Zt = 
1.55,1.03,0.66,0.46. 

As discussed in LC93 and LC94, the probability that a 
halo of mass M at redshift z has a progenitor with mass 
between M /2 and M at z p , gives the probability that its 
formation epoch was earlier than z p . In differential form, the 
probability distribution of formation redshifts is therefore 
given by: 



-?-(z c \M ,z ) = 
dzf 



M /2 



Mo 

M D 



9 ( d/ 



dzt VdM p 



dM n 



(7) 
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Figure 5. Differential probability distribution of formation red- 
shifts Zf for halos at z = with masses M /Mq = [10 8 , 5x10'], 
i = 11,12,13,14 (see labels), as a function of (1 + Zf). His- 
tograms and error bars result from the average over 10 realisa- 
tions of the MCM. N labels the number of halos found at z = 
in the corresponding mass range. Solid curves refer to the an- 
alytic prediction of the EPS theory computed from eq.(7), for 
M = Mo ~ 2xlO n , 2xl0 12 , 1.7xl0 13 , 1.4xl0 14 M Q . 



where d//dM p = d/(M p , Zf|M ,z )/dMp is the progenitor 
mass distribution (see Sec. 3.2). Formation times computed 
by means of the previous formula are found to be in good 
agreement with N-body simulation results, except for halos 
on cluster scales that form earlier than predicted by the 
EPS theory (LC94; Tormen 1998; note however that these 
conclusions are drawn for scale-free power spectra only). 

Once again, numerical results are derived for ranges 
of masses M , thus the analytic prediction for zt has been 
computed in the same way as discussed in Sec. 3.2. In this 
case however, the probability distribution for the average 
mass dp(z{\M z ) / dzf , and the average probability distri- 
bution (dp/dzt) are almost indistinguishable in the cho- 
sen range of M . Only the former is therefore shown in 
Fig. pi where solid curves correspond to eq.(7) solved for 
M = Mo ~ 2.2 x 10 11 , 2 x 10 12 , 1.7 x 10 13 , 1.4 x 10 14 M Q . 
Within the error bars, a very good agreement is found for 
the most massive objects. For intermediate mass halos, the 
epoch when they first appear, as well as the rising of the 
probability distribution with decreasing z is well reproduced 
by the MCM. However, they do not present the expected 
peak of formation epoch, but instead still form at very recent 
times, in contradiction with the expectations of the EPS the- 
ory. A severe disagreement is found for small objects, with 
the departure of MCM relative to the EPS theory going in 
the opposite sense. Low mass halos in fact preferentially col- 
lapse and stop forming at earlier epochs than predicted, with 
a peak of formation at about zt = 1.6, instead of zt — 0.85. 
No significant improvements are obtained if different values 
for the collapse threshold S c are adopted. This confirms once 
more that the history of low mass objects is not well followed 
in the model. 



3.5 Two— point Correlation Function 

Since the relative positions of halos within the box are known 
by construction, the MCM also contains information about 
their spatial distribution. We have computed the two-point 
autocorrelation function of DM halos, by counting the num- 



ber of objects separated by a distance r, and comparing it 
to the value expected for a Poissonian distribution: 



ar) 



Nnnjr) 

N RR (r) 



(8) 



where Ndd(j) is the number of pairs whose geometric cen- 
tres are separated by a distance between r and r + dr, and 
Nhr{t) is the same quantity if halos were randomly dis- 
tributed in the same volume: N RR (r) = (l/2)N 2 l (dV/V), 
where N Q is the total number of halos, dV is the volume of 
the shell at r with thickness dr, and V is the total volume 
of the box. 

Results for halos selected in four mass ranges at z = 0, 1 
and 3 are shown as circles in Figs. |6| Fig. Q , and Fig. |s|, re- 
spectively. Note that points are found at separations smaller 
than the typical halo sizes (marked by the vertical thick lines 
in plots) . This is a consequence of the nonspherical shape of 
halos in the MCM, allowing the distance between two cen- 
tres to be smaller than the spherical radius R artificially 
attributed to each object in this computation. 

Using an approach based on the EPS theory, CLMP 
give an analytic formula for the halo two-point correlation 
function, which is valid for separations r larger than R. In 
particular, when r 3>R, the correlation function of objects 
of mass M identified at redshift z, can be expressed as: 

Z hh (r,M,z)=b 2 1 (M,z)Z m (r,z) + h 2 2 (M,z)e m (r,z) + ...,(9) 

where £ m {r, z) is the matter correlation function (the Fourier 
transform of the power spectrum) linearly extrapolated to 
redshift z. The linear bias function bi(M, z) was already ob- 
tained by MW96 using a different approach still based on 
the EPS theory. It is given by: 

bi(M,z)= J° T -T 1 : ( 10 ) 



2 (M, 



where <r(M, z) is the mass variance linearly extrapolated to 
redshift z: a(M, z) - cr(M) (1 + z)" 1 . CLMP show that the 
second order bias factor is: 



&2(M,*) = 



1 



\M,z) 



2 (M, 



(11) 



If the typical nonlinear mass M,(z) for dark matter halos 
is defined as <7[M,(z), z] = S c , it results from eq.(10) that 
the first order bias vanishes for M=M*, and £hh is then 
determined by the second order term only. For redshifts 
z = 0, 1, 3, the values of M. are 3.4 x 10 13 , 1.2 x 10 12 , and 
6.2xl0 9 M Q , respectively. 

When a finite range of halo masses is considered in- 
stead of a single value of M, the theoretical halo correlation 
functions can still be estimated by eq.(9), with the two bias 
factors replaced by their mean values in the mass interval, 
weighted by the mass function n(M, z) = dn/dM: 



/ b t (M, z) n(M, z) dM 

Mi 

M a 

/ n(M,z)dM 



1,2. 



(12) 



Long dashed curves in Figs, kiL Fig. [j], and Fig. |s| have 
been computed by means of eq.(9)-(12) for the correspond- 
ing mass ranges and redshifts. Also shown for comparison are 
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the linear mass correlation function at each redshift (dotted 
lines), and £hh computed with the linear term in eq.(9) only 
(dashed-dotted lines). For all redshifts and halo masses, the 
autocorrelation function derived from the MCM is in a re- 
markable good agreement with the predictions of the EPS 
theory. Even if the analytic formula for £hh has been ob- 
tained in the limit of separations much larger than R, a rea- 
sonable agreement is also found when this condition is not 
exactly satisfied. Moreover, results of the MCM are well de- 
scribed by the linear bias relation also in the (slightly) non- 
linear clustering regime (i.e., for separations where £ m (>") 
is slightly larger than unity), and thanks to the second or- 
der term, eq.(9) still provides a very good description of the 
model correlation function, even for masses near M», where 
fei vanishes. 

Such an agreement between the MCM results and ana- 
lytic predictions derived from the EPS theory only ensures 
the reliability of the model in correctly taking into account 
the clustering of high mass (M £ M») halos, but it also high- 
lights its limitations for small objects. Indeed, accurate com- 
parisons with N-body simulations show that the correlation 
function given by eqs.(9)-(12) correctly describes numerical 
results for halos with masses larger than M„, but signif- 
icantly overestimates the clustering of small mass objects 
(PCL; Jing 1999; Sheth & Tormen 1999; Sheth, Mo, & Tor- 
men 1999). For M < M*, the analytic bias factor b\ is signif- 
icantly lower (more negative) than that found in numerical 
simulations, whereas a better fit to the N-body Lagrangian 
correlation function (with errors within the 15% for a CDM 
cosmology) is obtained by means of the linear term in eq.(9), 
with b\ replaced by (Jing 1999; but see also Sheth, Mo, & 
Tormen 1999): 



MM,*) 



q 4 (M, 
2 St 



+ 1 



(0.06-0.02n) 
X [l + b!(M,z)]-l, 



(13) 



where n is the index of the power spectrum P(k), computed 

as: 

_ dlnP(fe) 



dink 



\k=2-n/R- 



(14) 



The correlation functions computed as £;,h = bj £ m are 
plotted in Figs. pl-H as solid curves. Relative to these N- 
body based correlation functions, the MCM overestimates 
the clustering of halos on the low mass (M < M,) regime. 



Figure 6. Auto-corrclation function of halos with mass M se- 
lected at redshift z = 0. Each panel refers to a different mass 
range: M/Mq = [10* , 5x10"], with i = 11, 12, 13, 14 (see labels). 
Average results and standard deviation from 10 MCM realisations 
are plotted as circles and error bars, and the average number of 
halos found for each mass range is also indicated as N. Separa- 
tions are in units of the box length (L=100 Mpc = 256 bcs). The 
three vertical thick lines mark the typical Lagrangian radius R of 
halos in the given range of M: the two shorter ones correspond 
to the minimum and the maximum mass in the range, while the 
longer one refers to the mean mass in the interval, weighted by 
the mass function. Long dashed curves have been computed by 
use of the linear and the second order bias factors [Catelan et 
al. 1998; see eqs.(9)— (12)], and the values of bi and 62 are la- 
belled in each panel. Also shown are the linear mass correlation 
function at the given redshift (dotted curves), and ^^^ computed 
with the linear bias only, as first discussed by Mo & White (1996; 
dashed-dotted lines). Solid curves corresponds to the correlation 
functions computed as £/,/, = b j £ m , where the value of bj (see 
label) is derived from Jing's formula [eq. (13)], which provides a 
good fit to N-body simulations. The first order bias vanishes at 
M* ~3.4xl0 13 M Q . 



4 DISCUSSION 

As far as the halo mass function, and the conditional prob- 
ability distribution of progenitor and child halos are con- 
cerned, a good general agreement between MCM results and 
PS and EPS analytic formulae is found, but an underpro- 
duction of low- mass objects in the model is apparent. This 
limits the mass resolution of the MCM to a minimum of 
8 base cells (3 x 10 10 Mq , for the present choice of cosmo- 
logical parameters and box size). Compared to the mass 
function from N-body simulations (well described when a 
nonspherical model for the collapse of density fluctuations 
is considered; see LS98; LS99; Sheth, Mo, & Tormen 1999, 
and references therein), the MCM produces a significantly 
lower number of high-mass halos, especially at early times. 
Since a finite box is used for representing the Universe, the 



Figure 7. The same as in Fig. pi but for halos selected at z = 1. 
No results for halos with mass in the range [10 14 , 5xl0 14 ] Mq are 
shown because only 5 of them have already formed at this epoch. 
At z = l, M* ~ 1.2xlO 12 M . 



effective amplitude of the mass variance on large scales is 
smaller than that expected from the input power spectrum, 
used in the computation of the analytic formulae. Such an 
effect may in part be responsible for the underproduction 
of high-mass halos in the MCM with respect to theoretical 
predictions. Also changing the linking and the overlapping 
conditions (see Sec. 2.2) helps obtaining larger mass halos, 
especially at high z, but it is not clear which are the best 
criteria (see also Nagashima & Gouda 1997). Moreover, os- 
cillations in the mass functions occurring at block masses 
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Figure 8. The same as in Fig. H, but for halos selected at z = 
3. No results for larger masses are shown because too few, or 
no high-mass halos, have already formed at this redshift. Here, 
M» ~6.2xl0 9 M Q . 



of 8 1 base cells, i — 1,2,.., are apparent, but may possi- 
bly disappear if different criteria are adopted when deciding 
whether to merge or not pre-existing halos and form a new 
structure, as well as if set of grids displaced in a different 
way (no longer by half a block-length) are used. 

The distribution of formation redshifts is in good agree- 
ment with analytic predictions for high-mass halos, even if 
the peak of formation is systematically shifted towards more 
recent epochs. For intermediate mass objects, the MCM cor- 
rectly reproduces the analytic expectations only at high red- 
shifts, then it keeps forming halos also at very recent times, 
in contradiction with the EPS theory. Once more, a failure 
of the model in describing the history of low mass halos 
is evident, since they systematically form at earlier epochs 
than predicted. 

Finally, a remarkable agreement of the two-point cor- 
relation function is found with respect to the predictions 
derived from the EPS theory (MW96; CLMP), for all con- 
sidered masses and redshifts. This ensures that the model 
reliably retains information about the spatial correlation of 
high mass halos (M £ M»), but it also overestimates the 

clustering of small objects (as do all analytic formulae). The 
amplitude of their correlation function in fact is significantly 
(2-3 times) higher than that predicted by the fitting formula 
recently proposed by Jing (1999), that correctly describes 
the correlation function found in numerical simulations. As 
discussed by Jing (1998, 1999) and PCL, this difference be- 
tween N-body and EPS results in Lagrangian space, sug- 
gests that the criteria adopted in the PS theory for identi- 
fying bound virialised objects in the initial conditions are 
inadequate. The assumption of spherical symmetry for the 
collapse is certainly a strong simplification, and it also af- 
fects both the mass function and the typical formation epoch 
of structures. Considering that halos in the MCM are pro- 
duced with a large variety of shapes, we intend to adopt a 
nonspherical condition and study its effects on the resulting 
mass function and formation redshift distribution in a future 
paper. 

The physical processes ruling gas cooling, dissipative 
collapse, star formation, evolution and feedback (as well as 
interactions and merging between galaxies) are currently im- 
plemented in merging history trees of DM halos, so far ob- 
tained through two main approaches. In semi-analytic mod- 
els (KW93; Kauffmann, White & Guiderdoni 1993; Kauff- 



mann, Guiderdoni & White 1994; Cole et al. 1994; Baugh 
et al. 1998; Somerville & Primack 1999; and papers in these 
series), the merging history of DM halos is built through 
Monte-Carlo realisations of the block model or EPS formal- 
ism, with no or not accurate spatial information. More re- 
cently (Kauffmann, Nusser & Steinmetz 1997; Governato et 
al. 1998; Benson et al. 1999), DM halos have been selected 
from cosmological N-body simulations, but their merging 
trees are still computed with the Monte Carlo technique. In 
a "fully" hybrid model (Roukema et al. 1997; Kauffmann 
et al. 1999), merging trees are also computed from the out- 
put of large N-body simulations, and as a consequence they 
retain the spatial and dynamical information of the parent 
simulation, but they suffer from its limited mass resolution 
and expensive CPU cost. 

The interest of the MCM is that it represents an inter- 
mediate approach. It is very fast, and it partly retains spatial 
information in the linear or weakly non-linear regime. A pri- 
ori, it suffers from the same resolution problem as merging 
trees built from N-body simulations. For the same choice 
of cosmological parameters and box length, the 256 3 base 
cells have the same mass as the 256 3 particles, and reliable 
halos cannot be obtained below ~ 8 base cells or 10 par- 
ticles. However, its low cost in terms of CPU time allows 
to run realisations of sub-boxes, thus improving the mass 
resolution. Moreover, many choices of the cosmological pa- 
rameters, shape and normalisation of the power spectrum 
of linear fluctuations can be tested. So the MCM appears 
as a versatile and rapid method to test physical ideas about 
galaxy, group and cluster formation in various cosmologies, 
mostly when some degree of spatial information can be use- 
ful. 

In particular, the MCM can be suitable for studying 
galaxy clusters, mainly at low redshifts, where a good agree- 
ment between MCM and analytic results is found, not only 
in terms of mass functions, but also in the distribution of 
formation redshifts, as well as in the halo two-point corre- 
lation function. Also the population of Lyman-break galax- 
ies at z — 3 can be reasonably well studied by means of 
the MCM. In fact, these objects are often interpreted as 
star-forming galaxies located at the centre of halos of about 
10 12 M Q (e.g., Steidel et al. 1996; Giavalisco, Steidel & Mac- 
chetto 1996; Steidel et al. 1998; Giavalisco et al. 1998; Baugh 
et al. 1998; but see also Somerville, Primack & Faber 1998). 
For these masses and redshifts, the model provides a reason- 
ably good description of both the mass distribution and the 
formation history. Moreover, the correlation function fairly 
matches the numerical results over a large range of halo sep- 
arations, thus allowing in principle to investigate the clus- 
tering properties of Lyman-break galaxies. 



5 CONCLUSIONS 

The Merging Cell Model originally proposed by Rodrigues 

6 Thomas (1996) for a scale-free power spectrum, has been 
developed in the case of the SCDM cosmology. Its reliability 
has been tested not only in terms of the halo mass function, 
but also comparing the distributions of the progenitor and 
child masses, as well as that of halo formation times, to the 
analytic predictions derived by the Press & Schechter theory 
and its extensions. 



© 1994 RAS, MNRAS 000, [J]-|h] 



10 B. Lanzoni 



For the first time in the case of a semi-analytic merging 
tree model, we have also computed the halo two-point corre- 
lation function, and compared it to the available theoretical 
predictions. 

We have stressed the major successes of the model, as 
well as its main weakness, and several possible solutions to 
improve it have been proposed. 

Two main fields where the use of this method can be 
of particular interest have been recognised. It appears to be 
a suitable tool for studying the properties of cluster-scale 
objects, mainly at low redshift, as well as the population of 
Lyman-break galaxies, and their clustering at high z. 

We intend to apply the method in a more realistic cos- 
mological scenario (as the open and the lambda CDM), and 
directly test it against N-body simulations in a forthcoming 
work. 
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